Model Selection

Week 8

Evie Zhang

Old Dominion University

Topics

  • Seasonality
  • Model Selection
  • Exam Review

Rewind: AR(1)

\[Y_t = \alpha + \beta Y_{t-1} + e_t\]

\[E[Y_t] = \frac{\alpha}{1-\beta}\]

\[\text{Var}(Y_t) = \frac{\sigma^2}{1-\beta^2}\]

Rewind: AR(1)

\[Y_t = 10 + .3Y_{t-1} + e_t\]

Code 
shocks <- rnorm(200)
y <- c(shocks[1], rep(0, 199))
for(i in 2:200) y[i] <- 10 + .3*y[i-1] + shocks[i]

plot(y, type = "l", xlab = "", ylab = "",
     main = paste0("Mean: ", round(mean(y[100:200]), 1), "; Var: ", round(var(y[100:200]), 1)))

Rewind: AR(1)

\[Y_t = 10 + .3Y_{t-1} + e_t\]

Code 
shocks <- rnorm(200)
y <- c(shocks[1], rep(0, 199))
for(i in 2:200) y[i] <- 10 + .3*y[i-1] + shocks[i]

plot(y, type = "l", xlab = "", ylab = "",
     main = paste0("Mean: ", round(mean(y[100:200]), 1), "; Var: ", round(var(y[100:200]), 1)))

Rewind: AR(1)

\[Y_t = 10 + .3Y_{t-1} + e_t\]

Code 
shocks <- rnorm(200)
y <- c(shocks[1], rep(0, 199))
for(i in 2:200) y[i] <- 10 + .3*y[i-1] + shocks[i]

plot(y, type = "l", xlab = "", ylab = "",
     main = paste0("Mean: ", round(mean(y[100:200]), 1), "; Var: ", round(var(y[100:200]), 1)))

Rewind: AR(1)

\[Y_t = 10 + .3Y_{t-1} + e_t\]

Code 
auto.arima(y[100:200])
Series: y[100:200] 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1     mean
      0.4191  14.0641
s.e.  0.0907   0.1681

sigma^2 = 0.9958:  log likelihood = -142.18
AIC=290.37   AICc=290.62   BIC=298.21

Seasonality

Types of Seasonality:

  • Deterministic
    • arima(df$y, xreg = mat[,2:12])
  • Stochastic
    • arima(df$y, seasonal = list(order = c(1, 0, 0), period = 12))
    • arima(df$y, seasonal = list(order = c(0, 1, 0), period = 12))

Seasonality

How to get auto.arima to use seasonal components:

  • Suppose you have df$y, that is a monthly time series starting in 2005-06-01.
  • y_ts <- ts(df$y, start = c(2005, 06, 1), frequency = 12)
  • auto.arima(y_ts)

Seasonality

Code 
N <- 12*12
t <- 1:(N)
shocks <- rnorm(N)
v1 <- rep(0, N)
v2 <- rep(0, N)
v3 <- rep(0, N)

v1 <- t + shocks
v2[1] <- shocks[1] + 1
for(i in 2:length(t)) v2[i] <- .4*v2[i-1] + shocks[i] + 1
v3[1] <- shocks[1] + 1
for(i in 2:length(t)) v3[i] <- v3[i-1] + shocks[i] + 1

plot(t, v1, col = "tomato", type = "l",
     ylim = range(v1, v2, v3), xlab = "", ylab = "")
lines(t, v2, col = "dodgerblue")
lines(t, v3, col = "mediumseagreen")
legend("topleft", bty = "n", pch = 15, cex = 1.33,
       col = c("tomato", "dodgerblue", "mediumseagreen"),
       legend = c(expression(paste("Y"[t], " = t + e"[t])),
                  expression(paste("Y"[t], " = 1 + .4*Y"[t-1]," + e"[t])),
                  expression(paste("Y"[t], " = 1 + Y"[t-1]," + e"[t]))))
legend("left", bty = "n", pch = c(NA, rep(15, 3)), cex = 1.33,
       col = c("black", "tomato", "dodgerblue", "mediumseagreen"),
       legend = c(expression(paste("E[Y"[t=144],"]:")),
                  144, 1.67, 144))

Seasonality

\[Y_t = t + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v1); pacf(v1)

Seasonality

\[Y_t = 1 + .4Y_{t-1} + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v2); pacf(v2)

Seasonality

\[Y_t = 1 + Y_{t-1} + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v3); pacf(v3)

Seasonality

Code 
adf.test(v1)
adf.test(v3)

    Augmented Dickey-Fuller Test

data:  v1
Dickey-Fuller = -5.4424, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary


    Augmented Dickey-Fuller Test

data:  v3
Dickey-Fuller = -1.9447, Lag order = 5, p-value = 0.5995
alternative hypothesis: stationary

Seasonality

Code 
t <- rep(1:12, 12)
shocks <- rnorm(N)
v1 <- rep(0, N)
v2 <- rep(0, N)
v3 <- rep(0, N)

v1 <- abs(t-6) + shocks - 3
v12 <- abs(t-6) + shocks - 3 + (1:length(t)) / 12
v2[1:12] <- shocks[1:12]
for(i in 13:length(t)) v2[i] <- .4*v2[i-12] + shocks[i] + 1
v3[1:12] <- shocks[1:12]
for(i in 13:length(t)) v3[i] <- v3[i-12] + shocks[i] + 1

plot(1:length(t), v1, lwd = 3,
     col = alpha("tomato", .4), type = "l",
     ylim = range(v1, v2, v3, v12),
     xlab = "", ylab = "")
lines(1:length(t), v2, lwd = 3,
      col = alpha("dodgerblue", .4))
lines(1:length(t), v3, lwd = 3,
      col = alpha("mediumseagreen", .4))
lines(1:length(t), v12, lwd = 3,
      col = alpha("orchid", .4))
legend("topleft", bty = "n", pch = 15, cex = 1.33,
       col = c("tomato", "orchid", "dodgerblue", "mediumseagreen"),
       legend = c(expression(paste("Y"[t], " = ", alpha[S], " + e"[t])),
                  expression(paste("Y"[t], " = t/12 + ", alpha[S], " + e"[t])),
                  expression(paste("Y"[t], " = .4*Y"[t-12]," + e"[t])),
                  expression(paste("Y"[t], " = Y"[t-12]," + e"[t]))))

Seasonality

\[Y_t = \alpha_S + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v1)
pacf(v1)

Seasonality

\[Y_t = t/12 + \alpha_S + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v12)
pacf(v12)

Seasonality

\[Y_t = 1 + .4 * Y_{t-12} + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v2)
pacf(v2)

Seasonality

\[Y_t = 1 + Y_{t-12} + e_t\]

Code 
par(mfrow = c(1, 2))
acf(v3)
pacf(v3)

Seasonality

\[\Delta Y_t = Y_t - Y_{t-1} = 1 + (Y_{t-12} - Y_{t-13}) + e_t\]

Code 
par(mfrow = c(1, 2))
plot(diff(v3), type = "l")
acf(diff(v3))

Seasonality

\[\Delta_{12} Y_t = Y_t - Y_{t-12} = 1 + (Y_{t-12} - Y_{t-24}) + e_t\]

Code 
par(mfrow = c(1, 2))
plot(diff(v3, lag = 12), type = "l")
acf(diff(v3, lag = 12))

Seasonality

Code 
par(mfrow = c(1, 2))
plot(diff(diff(v3, lag = 12)), type = "l")
acf(diff(diff(v3, lag = 12)))

Seasonality Diagnostics

Is this stationary?

Look at the variance:

  • If increasing, stochastic (diff across period)

  • If constant:

    • If trend, deterministic (trend + FE)
    • If ACF/PACF is cyclic, deterministic (FE)
    • If ACF/PACF spikes, stochastic (AR(period))

Seasonality Diagnostics

\[Y_{t} = 1 + .4Y_{t-1} + .4Y_{t-12} + e_t\]

Code 
N <- 144*3
t <- rep(1:12, N/12)
shocks <- rnorm(N)
v4 <- rep(0, N)

v4[1] <- shocks[1]
for(i in 2:12) v4[i] <- .4*v4[i-1] + shocks[i] + 1
for(i in 13:length(t)) v4[i] <- .4*v4[i-1] + .4*v4[i-12] + shocks[i] + 1
plot(v4, type = "l")

Seasonality Diagnostics

Code 
par(mfrow = c(1, 2))
acf(v4); pacf(v4)

Seasonality Diagnostics

Code 
par(mfrow = c(1, 2))
acf(diff(v4)); pacf(diff(v4))

Seasonality Diagnostics

\[Y_{t} = 1 - .4Y_{t-1} + .4Y_{t-12} + e_t\]

Code 
v5 <- rep(0, N)

v5[1] <- shocks[1]
for(i in 2:12) v5[i] <- -.4*v5[i-1] + shocks[i] + 1
for(i in 13:length(t)) v5[i] <- -.4*v5[i-1] + .4*v5[i-12] + shocks[i] + 1
plot(v5, type = "l")

Seasonality Diagnostics

Code 
par(mfrow = c(1, 2))
acf(v5); pacf(v5)

Seasonality Diagnostics

\[Y_{t} = 1 + .4Y_{t-1} + .4Y_{t-12} + e_t\]

Code 
arima(v4, c(1, 0, 0),
      seasonal = list(order = c(1, 0, 0),
                      period = 12))

Call:
arima(x = v4, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))

Coefficients:
         ar1    sar1  intercept
      0.5013  0.5055     4.8013
s.e.  0.0455  0.0470     0.2038

sigma^2 estimated as 1.148:  log likelihood = -644.68,  aic = 1297.37

Seasonality Diagnostics

\[Y_{t} = 1 - .4Y_{t-1} + .4Y_{t-12} + e_t\]

Code 
arima(v5, c(1, 0, 0),
      seasonal = list(order = c(1, 0, 0),
                      period = 12))

Call:
arima(x = v5, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))

Coefficients:
          ar1    sar1  intercept
      -0.4558  0.4736     1.0223
s.e.   0.0450  0.0448     0.0625

sigma^2 estimated as 1.039:  log likelihood = -622.95,  aic = 1253.91

Forecasting Airplane Passengers

Code 
data("AirPassengers")
air <- AirPassengers
air <- log(air)
plot(air, ylab = "log(Passengers)")
abline(v = 120, lty = 2, col = alpha("black", .6))

Forecasting Airplane Passengers

Code 
plot(diff(air), ylab = expression(paste(Delta, "log(Passengers)")))

Forecasting Airplane Passengers

Code 
plot(diff(air, lag = 12), ylab = expression(paste(Delta[12], "log(Passengers)")))

Forecasting Airplane Passengers

Code 
plot(diff(diff(air, lag = 12)), ylab = expression(paste(Delta, " ", Delta[12], "log(Passengers)")))

Forecasting Airplane Passengers

Code 
par(mfrow = c(1, 2))
acf(diff(diff(air, lag = 12)))
pacf(diff(diff(air, lag = 12)))

Forecasting Airplane Passengers

Code 
air1 <- ts(air[1 : (length(air) - 24)],
           start = c(1949, 1, 0), frequency = 12)
reg <- auto.arima(air1)
reg
Series: air1 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.3424  -0.5405
s.e.   0.1009   0.0877

sigma^2 = 0.001432:  log likelihood = 197.51
AIC=-389.02   AICc=-388.78   BIC=-381

Forecasting Airplane Passengers

Code 
par(mfrow = c(1, 2))
acf(reg$residuals)
pacf(reg$residuals)

Forecasting Airplane Passengers

Code 
p <- predict(reg, n.ahead = 24)
plot(1:(length(air)),
     rep(0, length(air)),
     type = "n", xlab = "Time",
     ylab = "Air Passengers",
     ylim = exp(range(p$pred + 1.645*p$se, air1)))
lines(1:length(air), exp(air), col = "tomato", lwd = 2)
lines(1:length(air1), exp(air1), lwd = 2)
polygon(x = c(1:24 + length(air1),
              rev(1:24 + length(air1))),
        y = c(exp(p$pred + 1.645*p$se),
              rev(exp(p$pred - 1.645*p$se))),
        col = alpha("dodgerblue", .3), border = NA)
lines(1:24 + length(air1), exp(p$pred),
      col = "dodgerblue", lwd = 2, lty = 2)
abline(v = 120, lty = 2, col = alpha("black", .6))

Model Selection

Model Selection

  • “Economic” Theory
  • Measures of (in-sample) Fit
  • Repeated \(t\)/\(F\) tests
  • Bayes Information Criterion
    • selects the model most likely to be true (given the data)
  • Akaike Information Criterion
  • k-Fold Validation

BIC vs AIC

For both BIC and AIC, smaller is better.

Both BIC and AIC penalize the model when lags increase.

AIC penalizes less than BIC for more parameters, and therefore chooses models with more lags.

BIC picks which model is most likely to be true given the data whereas AIC picks\(^*\) the model with lowest forecast risk.

k-Fold Validation

k-Fold Validation

  • Select a handful of model specifications
  • Suppose 10 years of data.
    • Fit using years 1-5, forecast year 6
    • Fit using years 1-6, forecast year 7
    • …
  • Select the model with smallest (weighted?) average RMSE/MAPE/etc.

Forecasting Hazards in San Diego

ECON 707/807: Econometrics II